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Synchronization is a universal phenomenon found in many non-equilibrium systems. Much recent interest 
in this area has overlapped with the study of complex networks, where a major focus is determining how a 
system’s connectivity patterns affect the types of behavior that it can produce. Thus far, modeling efforts 
have focused on the tendency of networks of oscillators to mutually synchronize themselves, with less em¬ 
phasis on the effects of external driving. In this work we discuss the interplay between mutual and driven 
synchronization in networks of phase oscillators of the Kuramoto type, and explore how the structure and 
emergence of such states depends on the underlying network topology for simple random networks with a 
given degree distribution. We find a variety of interesting dynamical behaviors, including bifurcations and 
bistability patterns that are qualitatively different for heterogeneous and homogeneous networks, and which 
are separated by a Takens-Bogdanov-Cusp singularity in the parameter region where the coupling strength 
between oscillators is weak. Our analysis is connected to the underlying dynamics of oscillator clusters for 
important states and transitions. 


Collective behavior of complex networks is a very 
active field of theoretical and practical research. 
In particular, models of oscillator networks have 
drawn much attention due to their numerous ap¬ 
plications across diverse fields, with a particular 
emphasis on synchronization phenomena. Here 
we study the dynamics of coupled oscillators, 
subject to periodic forcing, on random networks 
with different degrees of connectivity, and un¬ 
cover many dynamical behaviors as a few param¬ 
eters are varied. We find that the unfolding of 
synchronized states, and the possibility of bista¬ 
bility among them, differs for networks depending 
on how heterogeneous the degree of local connec¬ 
tivity is. This is explained through a combination 
of analytic and numerical results. 

I. INTRODUCTION 

The tendency for populations of oscillators to synchro¬ 
nize their dynamics and produce large-scale collective os¬ 
cillations is relevant in a wide range of contextJIlf^. A 
particularly simple class of models for this behavior was 
proposed by Kuramoto, where each oscillator in a net¬ 
work is described by a phase variable, which has a ten¬ 
dency to oscillate at its natural frequency and in phase 
with its neighborJ^. This model has given insights into 
the dynamics of many systems, from the synchronization 
of coupled chemical oscillators and Josephson junction 
arrays, to correlations in visual cortex emeriments and 
coherence in neutrino flavor oscillations^HU^ 

Much recent work on the Kuramoto model has con¬ 
cerned synchronization on complex networks, where 
the transition to coherent oscillations depends on the 
properties of the network topolog^EH^. Some impor¬ 
tant results are vanishing synchronization thresholds 
and explosi ve tran sitions for networks with large degree 
flnctuation j^^l^^ l^. However, the effects of external driv¬ 
ing are much less known, and questions about how dif¬ 


ferent networks of oscillators respond to driving, and to 
what extent they can be controlled, have not been an¬ 
swered, even though in many circumstances, external 
fields are present!^. An important example is the net¬ 
work of pacemaker cells, which play a role in determin¬ 
ing mammalian circadia n rhythms, and can be driven by 
light-dark cyclej^ ^^ l ^^ l 

In what follows, we discuss the interplay between mu¬ 
tual and driven synehronization in random networks of 
phase oscillators with a given degree distribution. In 
particular, we present key aspects of the stability di¬ 
agram for the driven Kuramoto model on these net¬ 
works, focusing on the appearance of a codimension-three 
Takens-Bogdanov-Cusp singularity in the parameter re¬ 
gion where the coupling strength between oscillators is 
weak. This bifurcation description is used to explain 
various pathways to driven and mutual synchronization 
in terms of synchronized oscillator clusters and network 
topology. 

II. MEAN-FIELD REDUCTION AND ANALYSIS 

Kuramoto showed that a system of limit cycle oscil¬ 
lators, each near their own Hopf bifurcation, with weak 
coupling to their neighbors and fast amplitude equilibra¬ 
tion, have the following simple equations of motion: 

nn 

-^ = Wi + J y] Aij sin(6»j - 6»i), (1) 
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where Oi is the phase of the ith oscillator, with natural 
frequency uJi, coupling strength J, and adjacency matrix 
for the interaction network Aij. Under generic circum¬ 
stances (e.g., when the natural frequencies are randomly 
assigned according to a symmetric and unimodal distri¬ 
bution without correlations to the topology), this system 
undergoes a critical transition from incoherence to mu¬ 
tual synchronization once the coupling strength exceeds 
a threshold, resulting in a fraction of the network oscil- 
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lating at the average of the nat ural f requencies, and with 
a stationary phase distributiorP^^I^. 

A simple extension of the Kuramoto model that in¬ 
cludes a periodic driving force is given 

= uji ^ J '^^Aij sm{0j — 6 >i) + — Oi), (2) 
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with external field strength E and frequency Q. With 
similar assumptions for the natural frequency distribu¬ 
tion, we expect each term to have the following effects on 
the dynamics: the randomness in the frequencies causes 
oscillators to have disperse phases with monotonic build¬ 
up in time, the coupling tends to align the phases of 
neighbors in proportion to the number of connections 
in a local environment (which will vary across the net¬ 
work), and the driving field tends to force oscillators to 
move at the driving frequency and away from its natu¬ 
ral frequency. The interactions among these tendencies, 
both cooperative and competitive, will depend on the 
magnitude of each term and the network topology, and 
therefore we expect an intricate dynamics with multiple 
behaviors and transitions^. 


A. Degree class dynamics 


To clarify the dynamics, we attempt to find a reduced 
description of ([^. For convenience, we study the phases 
in the co-moving frame of the driving, (pi = Oi — Qt: 

^ + sm{4>j - (pi) - Esm{(pi), (3) 

i 

and consider random networks with a given degree distri¬ 
bution, p/c, that specifies the fraction of oscillators with k 
neighbors. In particular, we will study the annealed limit 
of random networks explicitly, for which Aij = 
where ki and kj are drawn from for a network of size 
N with average degree (/c), but note that our results are 
in qualitative agreement with quenched models (such as 
the configuration model, Fig@. For annealed networks 
we find oscillator dynamics: 


— Jkilm 
dt 


^ N{k) 


-h EXm\e 


( 4 ) 

from which we can define the complex order parameter 


or the average interaction strength (both magnitude and 
phase) that an oscillator feels along an edge to its neigh¬ 
bors. 

We are interested in the thermodynamic limit, N 
oo, in which it is useful to consider the density of oscilla¬ 
tors with phase (p at time t, given degree k and frequency 


cj, p(0, t; k^uj). This probability density satisfies a conti¬ 
nuity relation: 


dp _ d 
dt dp 


■ ^-icf> ^icf> 

I id — -j- ——— (Jkz-\r E) — (Jkz-\- E)j p 

\ 2% 2% J 


(6) 


with 


= J t] UJ, k)e^*(kod(p, (7) 


where g{uj) is the natural frequency distribution and z 
is the complex conjugate oi z. In order to solve (H we 
expand p into its Fourier components: 


p{(p,t;u,k) = ^ 1 


- '^an{t;u!,k)e 

n=l 


incf) 


+ C.C. 


(8) 


and look for simple power-series solutions of the form, 
an{t; UJ, k) = d^{t] uj, k) - an ansatz which was proposed 
by Ott and Antonsen, and that i s applic able in a wide 
array of Kuramoto model variantIn this case it 
gives the dynamics for a{t;uj, k): 

^ = i [Jkz + E]+i{Lo-n)a-t [jkz + E]a^, (9) 

Cit Zi z 

which completely specifies the order parameter: 

z{t) = J k)duj. (10) 


In addition, the dimensionality of the system can be fur¬ 
ther reduced by performing the natural frequency inte¬ 
gral, for which we assume: 


5(w) = 


TT 


7 _ 

(w - Wo)^ + 7 ^ 


( 11 ) 


a Cauchy distribution with median ujq and scale 7 . 
Generically, a{uj,k,t) has no poles in the upper-half of 
the complex cj-plane, and therefore we perform contour 
integration of ( [1Q| ) closed in this region^, which reduces 
the integral to the residue at the pole cjo + iy: 

z{t) = E 


where 

^ = p^Jkz + S] -(l + iA)afc-i[j7fcf + £’]a|, (13) 

with the dimensionless time, r = yt, and normalized pa¬ 
rameters: £ = F'/y, J = J/y, and A = (^2 — u;o)/y- 
This is the fundamental equation for the thermody¬ 
namic limit of the forced Kuramoto model on annealed 
networks. The dynamics has been reduced to a descrip¬ 
tion of the average contribution to the order parameter 
by nodes of degree k, with the size of the state-space 
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equal to twice the number of degree classes. In the fol¬ 
lowing, we will focus on networks that have power-law 
degree distributions with finite cutoffs and Poisson dis¬ 
tributions, 

= P^= fci ’ 

E k'-^ 

k'=l 

respectively, though analytic results are given for arbi¬ 
trary distributions. We will refer to the former as simply 
“power-law” for brevity, though the degree cutoff, Kcut^ 
will be specified when pertinent. In general, the cutoff 
determines the dimensionality of the reduced system, and 
we find that its value is relevant for heterogeneous net¬ 
work behavior, where large degree nodes can contribute 
significantly to the dynamics. 


B. Limiting states 


First, we consider the states of mutual and driven syn¬ 
chronization in instructive limits. For instance, in the 
limit where f ^ 0, Eq.(p!^ describes an un-driven net¬ 
work, and has stable solutions corresponding to oscillat¬ 
ing waves: ak = r/c(T)e“*^^, z = 

^ = ^JkR[l - rl] - rk, (15) 

which reproduces known result^. In particular for the 
frame where = 0, the network tends to a purely os¬ 
cillating state at the average natural frequency cjo, with 
some fixed r^. In addition, the incoherent state, = 0, 
has a linear stability exponent 


2{k) 


(16) 


which implies a threshold for the onset of mutual synchro- 

nization in the absence of driving, ^ 0, 2{k) ~ 

consider situations where J > Jc^ and a cohere nt mu tu- 
ally synchronized state is stable without forcin^^^^. 

On the other hand, in the limit where the driving fre¬ 
quency is equal to the average natural frequency, A ^ 0, 
Eq.(13) describes states of driven synchronization, where 
the network is oscillating at the driving frequency on av¬ 
erage, with amplitudes given by the fixed points of the 
self-consistent equation: 


T.* -1 + \/l+ {JkR* +£f \ 

^ {k)\ {jkR*+£) y 


(17) 


C. Partial stability diagram 


Next, we provide the results of a stability and partiaP^ 
bifurcation analysis that delineates the boundaries be¬ 
tween the limiting states and helps to explain how each 
can be converted into the other. The associated stability 
diagrams are somewhat complicated, and it is therefore 
useful to have the results in hand before proceeding to 
fill in the details. A quantitative discussion and analysis 
can be found in Sec |IID| that derives some of the results, 
with a broader summary of behaviors found in Sec III A 



FIG. 1. (Color online) Stability diagrams for the driven Ku- 
ramoto model on random networks shown as functions of the 
driving field strength 8 and frequency detuning A. A legend 
is given in the bottom left, (a) Schematic diagram (left) for 
homogeneous behavior (e.g., power-laws with large exponents, 
k-regular, and Poisson degree distributions); (right) diagram 
for a power-law network with = 2, s = 3.0, and Kent = 200. 
(b) Schematic partiaP^ diagram (left) for heterogeneous be¬ 
havior (power-law degree distribution with sm all ex ponent), 
in the parameter region where J is weak (Sec jll D| ; partial 
diagram for a power-law network with J = 0.25, s = 2.3, and 
Kent — 1000, shown in two parameter ranges for clarity. A 
table indicating possible states is shown in the bottom right 
for important regions (Roman numerals). 


and with a large number of nodes entrained to the driv¬ 
ing. From 0 we can see that incoherence is not a so¬ 
lution when f 7 ^ 0, meaning that external driving al¬ 
ways enforces some level of coherent oscillations at its 
frequency. Moreover, multiple coherent solutions can ex¬ 
ist depending on the parameter values. 


The schematic stability diagrams shown in Figj^ illus¬ 
trate two types of behavio r in t he (A, 8) plane when the 
coupling, is weak (Sec |IID[ ). If we consider networks 
with power-law degree distributions, Fig§a) shows the 
generic behavior when the degree exponent, s, is large. 
We find that this is maintained for networks with rela- 
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tively homogeneous degree distributions, such as Erdos 
Renyi, k-regular, or complete graph^^. Conversely when 
the degree exponent is small, i.e., the degree distribu¬ 
tion has a heavy tail, the behavior looks like Figj^b). 
Because the former reproduces the behavior for the com¬ 
plete graph, and the latter occurs as the amount of varia¬ 
tion in the degree distribution is increased, we distinguish 
these cases by the terms homogeneous and heterogeneous 
driven behavior. 


In particular, we find that the Takens-Bogdanov point 
appears on the upper branch of the saddle-node bifur¬ 
cations for the homogeneous case, but appears on the 
lower branch for the heterogenous case, as depicted in 
FiglU The transition between the two behaviors can oc¬ 
cur, for example, by decreasing the degree exponent, for 
some fixed J and K^uti until the Takens-Bogdanov and 
cusp bifurcations are coinc ident , which typically occurs 
for some 2 < s < 3 (see Sec IID for descriptions of these 
bifurcations). The existence of this singularity allows us 
to construct the behaviors shown through a combination 
of analytic results, numerical continuation, and general 
predictions for subsequent bifurcations. Details are given 
in the following sections, and example stability diagrams 
are shown alongside schematics in Figj^ 


D. Stability analysis and bifurcations 


We begin constructing the stability digrams by first 
finding the fixed points of Eq.(13), which denote states 


of driven synchronization, and establish how such states 
change stability. In general, fixed points satisfy = 0, 
which implies the self-consistent condition for 


_ ^Pk 

~^W 


_ kpk “ (1 + + y^(l + + \Jkz* -h f p 


Jkz* £ 


(18) 


Each ak is a complex number, and so could be repre¬ 
sented by a magnitude and phase, or with real and imagi¬ 
nary parts. Next, it is useful to consider how the dynam¬ 
ics respond to perturbations away from the steady-states 
given by ( p^ , e.g. lZe[al\ -\-Xk and Xm[a^] -\-yk where 
and Uk are the /c’th components of the right eigenvectors 
of (13) at a^, in the real and imaginary part representa¬ 
tion of a/c. Equivalently, we can define r]k = ~^i^k -^Wk) 

and ffk = -^{^k ~Wk)^ with the perturbations + v^77/e 

and al -|- V2r]k . It is more convenient to use the latter and 
leave Eq.([T^ in its complex form, while keeping in mind 
that the standard results of bifurcation theory pertain to 
some underlying real representation of (13). 


We look for the linear stability spectrum of eigen- 
modes around a fixed point by adding the perturbations 


discussed into ( |13[ ), and collecting terms of order y: 

.,2 

k 


dr]k 

dr 


2 


E k'pk' ^2 k'pk' 


(k) 


QkVk^ 


dr]k 

dr 


with 


2 


E k'pk' ^ 

liT*' 


dh 


E k'pk' 


ql = l+iA + {Jkr+£)al; 


QkVk, 

(19) 

( 20 ) 


This system has a set of solutions, ^ = Xr]k and = 
Xrjk, from which we can find a self-consistent equation 
for the spectrum {A}. Solving for ijk and ffk in (19), 
multiplying by summing over /c, and eliminating the 


kpk ; 


constants -^r]k and Y,k "^dk gives: 


E 


>!Tk pkCLk 
2(/c)(A+q^) 


E 


Jk^Pk 

2{k){\+ql) 


E 


>!Tk pkCLk 
2(A:)(A+qJ) 


E 


Jk^Pk 

2{k){\+q*) 


= 1 - ( 21 ) 


Next, we catalogue relevant bifurcations found in 
Figpl and discuss their dynamical behaviors in Sec. 
ra First, the spectrum condition can be used to find 
the local codimension-one bifurcations, where some num¬ 
ber of eigenvalues cross the imaginary axis (codimen¬ 
sion implying the number of parameters th at mu st be 
changed in order for a bifurcation to occurThe 
most generic such crossing is the saddle-node bifurca¬ 
tion (SN), in which the spectrum at the equilibrium has 
one simple zero eigenvalue, and at which two equilibrium 
points collide and disappear: 



*4” k PkCLk 
‘^{k)ql 


Jk^Pk _ 1 
2{k)ql ^ 


2 


= 1 . 


( 22 ) 


The SN condition ( [E| ) predicts when steady states of 
driven synchronization vanish, and signifies when a lo¬ 
cal barrier (represented by the saddle) in the dynamics 
has been overcome. Importantly, we find that the lower 
branch of SN bifurcations contains a section of saddle- 
node-infinite-period bifurcations (SNIPER) (e.g., cross¬ 
ing V-III in Fig{^ where an SN occurs on a limit cycle 
of infinite periocP^. 

Another local codimension-one bifurcation is the Hopf 
(i^), in which the spectrum at the equilibrium has two 
purely imaginary eigenvalues, with all others having non¬ 
zero real parts. At this point the amplitude of a periodic 
orbit decreases continuously to zero with its period tend¬ 
ing to 27rI ujh^ where A = iujn'- 


2 

E Jk Pk<^k 
^ 2{k){iujH-\-ql) 

E Jk'^Pk _ 

^ 2{k){iujHrql) 


2 ^ 

E J'k Pk^k 
^ 2{k)iiujH-\-ql) 

E Jk^Pk _ 

^ 2{k)(iujH+ql) 


1. (23) 
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When the periodic orbit associated with the Hopf bifur¬ 
cation is stable, it is called supercritical (Hs uv)-, a nd when 
it is unstable, it is called subcritical {Hgub) In con¬ 

trast with homogenous network behavior (e.g., crossing 
la — Ilia in FigQa))!^, both branches of cycle-stability 
can appear if the degree distribution is broad enough 
(e.g., crossing II 5 — III 5 in Figj^b)). 

Beyond the local codimension-one bifurcations, there 
are two key local codimension-two bifurcations. These 
are important to unravel because they can inform us as 
to what global bifurcations occur. The first appears when 
two branches of the SN collide, in the neighborhood of 
which there exist three states of driven entrainment; this 
is known as a cusp (C) To find the C point, we first 
consider that near a bifurcation, the equations of mo¬ 
tion can be restricted to a center manifold with the same 
dimension as the number of eigenvectors whose eigen¬ 
values cross the imaginary axis, and is tangent to those 
vectors. Furthermore, the dynamics of the center mani¬ 
fold are equivalent to the normal form for the bifurcation. 
In the simple case of a SN^ the center manifold is one¬ 
dimensional, m = WT] + w‘^h + 0{w^) with the normal 
form: ^ = cw‘^ + 0{w^) 

The C bifurcation occurs when c = 0 , a condition for 
which can be found by substituting the center manifold 
expansion and normal form into (13), collecting terms 
of order and taking the complex inner product of 
the resulting vector, B, with the left eigenvectors, Q 
The right and left eigenvectors are found from a similar 
self-consistent analysis as for ( 21 ), and in the complex 
representation are respectively: 


AJk f x{X)-af 

mW - — I j + 


VkW = 


Cfe(^) = 


a (A) 


AJk i I — x{X)c 


(fc) 


A + ^ 


Zkpk ( _ x{\) _ 

<‘> V(E 4 Sfe- 0 (^ + «) 


_ Zkpk ( 




(24) 

(25) 

(26) 

(27) 


with constants A and Z, and with the conveniently de¬ 
fined sum. 


E 


x(A) = 


^k P 

2{k)(\+ql) 


E 


Jk^Pk 

2{k}{X+ql) 


- 1 


(28) 


Collecting terms of order in the expansion produces 
the bilinear form for (13) evaluated at the vector rik,rik- 


Bk{X) = -2JAkpkal - (Jkr + S)pl (29) 
Bk{X) = -2JAxk^kK - (Jkz* + £)7ii. (30) 


Putting these together generates the normal form coeffi¬ 
cient c, and a condition for the cusp bifurcation: 

c = 5ICfe(0)Bfe(0) + Cfe(0)Bfc(0) = 0, (31) 

k 

in conjunction with ( 20 ). 

It should be noted that for power-law networks with 
J < 2.5, the number of possible fixed points for this 
system is three, which we call the weak coupling re¬ 
gion. However, when the coupling is stronger, a degen¬ 
erate C point seems to emerge, which generates addi¬ 
tional unstable and saddle states, and complicates the 
unfolding (shown in FigQ, though much of the general 
structure is maintained for larger J. In this work, we re¬ 
strict ourselves to the weak coupling region for power-law 
networks, because the comparison between homogeneous 
and heterogeneous graphs is more straightforward. 

The second local codimension-two bifurcation is the 
Takens-Bogdanov (TH), at which the spectrum has a 
double root at zero. Attached to this bifurcation are 
curves of 5'A^ and H bifurcations as well as a curve of 
Homoclinic Bifurcations {HC) latter, the 

period of a cycle diverges as it collides with a saddle- 
point and connects its stable and unstable manifolds 
(e.g., crossing IV^ — in Figj^a)). To find the location 
of the TB bifurcation, we expand ( [2^ in powers of ujh^ 
and enforce that terms of order ujh vanish, which gives 
the criterion: 


IZe 


E 


Jk^Pkat 


E 


2{k)qf 

Jk'^Pk 


^{k)qt 



JkhkK 

2{k)rk 


xJk^Pk , 


= 0, (32) 


that in conjunction with ( 22 ), determines the bifurcation 
point. 

Finally, the highest codimension bifurcation that we 
consider arises when the C collides with the TB (TBC)^ 
implying that (22), (31), and (32) are all satis fied (wh ich 
also occurs in the Hodgkin-Huxley equationsIn 
addition to the bifurcations discussed, this particular sin¬ 
gularity predicts curves of codimension-two homoclinic 
bifurcations to Saddle-Node-Loops (SNL) and Neutral 
Saddles (NS), and curves of Degenerate Hopf bifurca¬ 
tions (DH). The latter two are termination points for 
curves of Limit-Point-of-Cycles {LPC). These bifurca¬ 
tions imply new behaviors that do not appear for ho¬ 
mogeneous networks and have interesting effects on the 
dynamics. Specifically, the LPC transition entails that 
a stable cycle collides with an unstable cycle and disap¬ 
pears, while the DH entails that an LPC emerges on 
a H point - typically as the Lyapunov exponent of the 
Hopf cycle vanishes. Because these bifurcations only oc¬ 
cur when the TH is on the lower branch of the SN^ they 
are not seen in homogeneous networks. Lastly, the SNL 
and NS entail that a homoclinic cycle is coincident with 
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di SN and a saddle whose whose eigenvalues sum to zero, 


III. OSCILLATOR DYNAMICS 


In this section, we explore some of the implications of 
the behaviors discussed on the dynamics of driven homo¬ 
geneous and heterogeneous networks. Both the mean- 
field (Eq|T^ and oscillator (EqJ^ dynamics are exam¬ 
ined. 


A. Key transitions and bistability 


Eirst, we can distill from the above that there are three 
primary ways that a stable mutually synchronized state 
can be created: Hgup^ SNIPER^ and LPC transitions. 
Qualitatively, we can think of such states as limit cy¬ 
cles, and can consider how their average amplitude and 
frequency (inverse period) emerge through each transi¬ 
tion. If we imagine changing one parameter (e.g. f), 
one of three things happens: the amplitude can appear 
continuously with a discontinuous frequency {Hgup)^ the 
amplitude can appear discontinuously with continuous 
frequency {SNIPER), or the amplitude and frequency 
can both appear discontinuously (LPC). The special 
case of continuous amplitude and frequency appearance 
occurs through a T5 bifurcation. Eigj^ shows a com¬ 
parison between the behaviors of mutually synchronized 
states produced by crossing these transitions. 

Interestingly, we find that each transition has a signa¬ 
ture in the average phase build-up with respect to the 
driving field. Eor example if the SNIPER transition 
is crossed (e.g., crossing V — III in EigQ, the order- 
parameter dynamics is a large limit cycle that includes 
the origirP^. This implies that the average phase of the 
network grows monotonically with respect to the field, 
and is therefore largely de-pinned from it, with a macro¬ 
scopic number of nodes lapping it continually. Moreover, 
this behavior holds widely for degree classes as well - 
most degrees continually lap the field on average, perhaps 
excluding low degree nodes (e.g., k=l or 2 ) depending 
on the parameters (EigJ^b)). On the other hand if the 
Hsup is crossed (e.g., crossing I — III in Eigffl, a small 
limit cycle emerges, centered around an unstable driven 
state. In this case there is no net build-up of the average 
phase with respect to the field; the motion is analogous 
to quasi-periodicity with average frequency equal to the 
driving, and an emergent “wobble” frequency given by 
(23) “ This behavior holds for all degree classes, im¬ 


plying that large and small degrees on average both have 
phase-trapped dynamics (Eigj^(a)). However if the LPC 
transition is crossed (e.g., crossing I 5 — II 5 in Eigj^b)), 
a large cycle emerges for the order-parameter that in¬ 
cludes the origin (similar to the SNIPER), but only holds 
for nodes with large degree on average, i.e., nodes of 
small degree undergo phase-trapped motion, while nodes 



FIG. 2. (Color online) Comparison of mutually synchronized 
states that arise from perturbations to driven states just be¬ 
low the key transitions for networks with power-law degree 
distributions. Subplots (a-c) show the average phase deflec¬ 
tion from a driven state for various degree classes versus time; 
colors for every degree class are specified in (b). (a) Below 
a Hsup-. S = Sh - 10“^, A = 5.5, J = 0.25, s = 2.6, and 
Ecut = 1000; the state appears with finite frequency and 
small amplitude, (b) Below a SNIPER: S = Ssn — 10“^, 
A = 0.2, H = 0.25, s = 2.3, and Kent = 1000; the state ap¬ 
pears with a large period and with all degree classes increas¬ 
ing phase monotonically with respect to the field, (c) Below 
a LPC (below Hsub)' S = Sh - 10-^ A = 5.5, J = 0.25, 
s = 2.3, and Kent = 1000; the state appears with finite fre¬ 
quency and large amplitude, and with high degree nodes in¬ 
creasing phase monotonically with respect to the field (phase- 
slip motion), while small degree nodes remain phase-trapped 
(on average), (d) Mutually synchronized state of with 
(c) parameter values, illustrating the cycle size variation with 
degree for states produced by crossing the LPC transition. 


of large degree undergo phase-slip motion (EigJ^c)). If 
we consider moving up the LPC by increasing £, more 
and more high degree nodes become trapped by the field, 
until all are trapped, and the Hsup occurs - the opposite 
limit brings us to the lower SNIPER (see FigQ- 

Another important difference between heterogeneous 
and homogeneous behavior concerns the bistability of 
driven and mutually synchronized states. Phase portraits 
are given in Fig 0 projected onto the order parameter, 
which demonstrate the behavior in important parame¬ 
ter regions. Eor homogeneous networks, bistability exists 
in a small region of parameter space, confined between 
the C and HC bifurcations (i.e., regions 11^ and IVa in 
EigQa)). In this case, there is bistability between two 
states of driven synchronization (region 11^ and EigJ^a)), 
until the Hsup is crossed (e.g., crossing 11 ^^ — IV^), and 
bistability between a state of quasi-periodic mutual syn¬ 
chronization and driven synchronization (Eigj^b))P^. In 
both cases, the manifolds of the saddle act as a separatrix 
between the two stable states. In contrast, for heteroge¬ 
neous networks there is only bistability between a large- 
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amplitude state of mutual syuchrouizatiou aud a siugle 
state of driveu syuchrouizatiou. The mutually syuchro- 
uized state eucloses all three fixed poiuts iu regiou IV 5 
(Figj^b) aud FigJ^d)), aud exists iu au additioual re¬ 
giou that does uot coutaiu a saddle ( II 5 iu Figj^b) aud 
FigJ^c)). Au example comparisou of the bistability iu 
fiuite uetwork simulatious for the two types of behavior 
is showu iu Fig|^ 


(a) 





C J 




FIG. 3. (Color online) Phase portraits of the dynamics pro¬ 
jected onto the complex plane of the order parameter for 
bistability regions (Ila, IVa, IF, and IVb, shown in (a),(b),(c), 
and (d), respectively). Initial transients were ignored, and 
curves were plotted once an effective two-dimensional dynam¬ 
ics was seen. The colors red and blue denote stable and un¬ 
stable fixed points and cycles, respectively. Panels have been 
rotated and scaled for clarity. 


B. Cluster behavior 

Finally, we are interested in how the states and transi¬ 
tions discussed in the previous sections appear at a finer 
scale of resolution: the dynamics of oscillator clusters 
in the network. For stable states of driven synchroniza¬ 
tion, we find a single macroscopic cluster of phase-locked 
nodes, which are entrained to the driving and are sta¬ 
tionary in the co-moving frame (labeled “L” in Figj^a)). 
This cluster is comprised of oscillators that have natural 
frequencies near the driving, with frequency ranges for 
degree classes that typically increase with degree, so that 
higher degree classes are able to stabilize a broader range 
of frequencies. Moreover, nodes with natural frequencies 
outside of their degree class’s locked range have average 
velocities (time average of Eq. (§) that are monotoni- 
cally increasing with the displacement from that range, 
and thus lap the driving field continually with disperse 
phases from one another. We therefore call these oscil¬ 
lators “winding” (labeled “W” in Figj^a)). The average 
velocities of phase-locked and winding nodes are shown 



FIG. 4. (Golor online) Gomparison of bistability of mu¬ 
tual and driven synchronization on heterogeneous and ho¬ 
mogeneous networks, (a) Stable states for a network with 
power-law degree distribution: ^ = 7.1, A = 8.0, J — 2.0, 
s — 2.0, and K^ut — 200; the order parameter is shown for 
the mean-field cycle (blue/black), annealed cycle (green/dark 
grey), configuration-model cycle (red/medium gray), mean- 
field equilibrium (magenta/medium gray point), annealed 
equilibrium (cyan/light gray), and configuration-model equi¬ 
librium (yellow with triangle/light gray), with good agree¬ 
ment among the respective states (region IF in Fig. [i](b)). 
Networks consist of 30,000 nodes, (b) Analogous plot for a 
Poisson degree distribution network with the same average 
degree as (a) and with J — 0.75; E — 2.27, A = 2.1216 for 
the mean-field, and 8 = 2.3, A = 2.134 for the annealed (re¬ 
gion IVa in Fig. [^(a))^. The arrow indicates where (b) can 
be found in z’s complex plane for comparison with (a). 


in Figj^a) for a driven state as functions of their degrees 
and natural frequencies. 

On the other hand, a stable state of mutual synchro¬ 
nization has a macroscopic cluster of nodes which lap 
the field together at some emergent average velocity. In 
addition, there exist other large “plateau” clusters of 
higher harmonics with average velocities that are inte¬ 
ger multiples of the fundamental velocity, and therefore 
lap the driving field 2,3,4... times in one network cycle 
(labeled as 1 , 2 ... in FigJ^b)). Collectively these har¬ 
monic plateaus drive phase-trapped nodes at a frequency 
equal to the fundamental velocity, causing them to wob¬ 
ble around the driving-field, but with average velocity 
zero (labeled 0 in Figj^b)). The last group of oscillators, 
which are between the plateaus, wind with average veloc¬ 
ities that grow monotonically with the displacement from 
a given plateau, and have disperse phases. This picture 
is consistent with general results for Kuramoto models, 
in which devil’s staircases do not appe ar, and velocities 
strictly increase between plateau^^^^^. Figgb) shows a 
typical velocity profile for a mutually synchronized state. 

Since we find that an important difference between 
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FIG. 5. (Color online) The average velocities for a network 
with a power-law degree distribution shown versus natural 
frequency, (a) A stable state of driven synchronization, (b) 
A stable state of mutual synchronization. The parameters are 
S = 7.1, A = 8.0, J = 2.0, s = 2.0, Kent = 200, and N = 
30000 at which (a) or (b) can be realized, given appropriate 
initial conditions (region Ilb in Figj^b)). The inset panels for 
(b) show 4)1 vs. T over one cycle with ranges [0,-7r], [0,-27r], 
[0,-47r], and [0,-67r] for the plateau numbers n — 0,1,2, and 
3, respectively. The node types for (a) are labeled “L” for 
locked “W” for winding, and appear next to the inset panels. 
Also, arrows indicate which cluster of oscillators are shown. 
A color legend for degree is given in (a). 


driven and mutually synchronized states is the appear¬ 
ance of plateaus in the velocity profile, we would like 
know how the plateaus are occupied when crossing the 
key transitions^. For instance, in crossing over a 
SNIPER transition, we find that the plateaus of the mu¬ 
tually synchronized state emerge from the phase-locked 
cluster of a driven state, as a finite fraction of locked 
nodes with natural frequencies near the average break 
away from the external field (shown in FigJ^a-b)). This 
is consistent with the de-pinning, continuous frequency 
and discontinuous amplitude appearance predicted by 
the mean-field dynamics. Conversely when crossing over 
the Hsup^ we find that a small stable cluster of wind¬ 
ing nodes in a driven state, with natural frequencies 
near the average, coalesce around the same average ve¬ 
locity ( [^ , and form the first plateau. As the transi¬ 
tion is approached the size of each plateau goes to zero 
(shown in FigJ^c-d)). This produces the discontinu¬ 
ous frequency and continuous amplitude limit cycle with 
quasi-periodicity described by the mean-field H bifurca¬ 
tion. Different still, when crossing over the LPC^ we find 
that a large group of nodes, which could form the wind¬ 
ing and locked clusters of a driven state, can coalesce 
around an average velocity instead (given appropriate 
initial conditions). This produces an additional stable 
state of mutual synchronization that is bistable with the 
driven state, and has plateaus that are disproportionately 



(C) 


0.1 

0.0 



FIG. 6. (Golor online) (a) A histogram of the average speeds 
for a network with a power-law degree distribution just below 
a SNIPER transition. The parameters are 8 = 3.61, A = 
4.0, J = 2.0, s = 2.0, Kent = 200, and N = 30000. We can 
see that roughly 30% of the network becomes de-pinned and 
forms the first plateau, which has speed (| (^) |n=i ^ 0). 
(b) The magnitude of the order parameter as a function of 
time, displaying relaxation dynamics as the network slows 
down in the neighborhood of a driven state that vanished in 
the transition, (c) Analogous histogram for a Poisson degree 
distribution network with the same average degree as (a), but 
just below a Hsup transition with parameters 8 = 2.41, A = 
2.6, H = 0.75. In this case, only 5% of the network occupies 
the first plateau, which has non-zero speed as the transition 
is approached, (d) Analogous plot to (b), showing the small 
amplitude, fast dynamics produced by crossing the Hsup- 



occupied by high degree classes. The velocities and order 
parameter dynamics are compared in Figj^and FigQa) 
for these bistable states, respectively. 


IV. CONCLUSION 


In this work we have studied the periodically driven 
Kuramoto model on random networks with a given de¬ 
gree distribution. A low-dimensional description was 
found, and a stability and partial bifurcation analysis de¬ 
veloped, which allowed us to predict many of the states 
and transitions of the model for sufficiently weak cou¬ 
pling between node^. In particular we found a Takens- 
Bogdanov-Cusp (TBC) singularity, appearing for power- 
law degree distribution networks as the degree exponent 
was lowered, which separated branches of heterogeneous 
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and homogeneous network behavior. The unfolding of 
this singularity was used to uncover important dynam¬ 
ical transitions including: Saddle-Node-Infinite-Period, 
Hopf, and Limit-Point-of-Cycles (LPC), as well as mul¬ 
tiple bistability regions that differed for the network 
types. Interestingly, we found that heterogeneous net¬ 
works do not support bistability of driven synchronized 
states or bistability of quasi-periodic synchronized states 
and driven states (which is the case for homogeneous 
networks), but only bistability of large amplitude mu¬ 
tually synchronized and driven states. Moreover, we dis¬ 
covered that the LPC transition for the heterogeneous 
branch occurs with phase-slip dynamics for nodes with 
high degree and phase-trapped dynamics for nodes with 
low degree (on average), implying a new route to mutual 
synchronization for driven heterogeneous networks which 
allows for qualitatively different behavior depending on a 
node’s degree. In addition, the structure of synchroniza¬ 
tion clusters for mutual and driven states was discussed 
and their transitions associated with bifurcations. 

Still, we have yet to resolve all of the transitions as¬ 
sociated with unstable cycles in the heterogeneous case 
(which could inform other interesting features of the dy¬ 
namics), and the full unfolding of network bifurcations 
in the strong coupling region. Moreover, many real net¬ 
works of interest have richer architecture than the sim¬ 
ple degree heterogeneity discussed he re: su ch as modu¬ 
lar, fractal, and multi-scale structnr j^^ l ^^ i The effects 
of these features on network synchronization are inter¬ 
esting subjects for future work. Finally, the control of 
complex networks is of immense interest, both theoreti¬ 
cal and practical. Our results can offer insight into the 
problem of controlling disordered oscillator networks. 
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